Characterization of lab-based swarms of Anopheles gambiae mosquitoes using 3D-video tracking

Mosquito copulation is a crucial determinant of its capacity to transmit malaria-causing Plasmodium parasites as well as underpinning several highly-anticipated vector control methodologies such as gene drive and sterile insect technique. For the anopheline mosquitoes responsible for African malaria transmission, mating takes place within crepuscular male swarms which females enter solely to mate. However, the mechanisms that regulate swarm structure or that govern mate choice remain opaque. We used 3D-video tracking approaches and computer vision algorithms developed for the study of other complex biological systems to document swarming behavior of a lab-adapted Anopheles gambiae line in a lab-based setting. By reconstructing trajectories of individual mosquitoes lasting up to 15.88 s, in swarms containing upwards of 200 participants, we documented swarm-like behavior in both males and females. In single sex swarms, encounters between individuals were fleeting (< 0.75 s). By contrast, in mixed swarms, we were able to detect 79 ‘brief encounters’ (> 0.75 s; < 2.5 s) and 17 longer-lived encounters (> 2.5 s). We also documented several examples of apparent male-male mating competition. These findings represent the first steps towards a more detailed and quantitative description of swarming and courtship behavior in one of the most important vectors of malaria.


A toy model: uniform circular motion.
In Section 2.2 of the main article, we provide a characterization of male mosquitoes' motion based on the analysis of the velocity probability distributions, and of the autocorrelation functions shown in Fig.4, Fig.5 and Fig.6: males essentially fly at a fixed altitude, performing consecutive, ring-like patterns in horizontal planes parallel to the ground. The narrow probability distribution of the y-component of the velocity, the one parallel to the direction of gravity shown in Fig.4a, indicates little deviation of v y from 0 or, equivalently, only slight variation of the direction of motion with or against the gravity, hence an almost constant flight altitude. On the other hand, the double-peaked distributions of v x and v z , shown at a group level in Fig.4b and Fig.4c and at an individual level in Fig.6 (third column), suggest a pseudo-circular motion in the horizontal plane. This correspondence between double peaked distribution and circular motion might not be intuitive to a general audience. Therefore, to aid understanding of the results presented in Section 2.2, we introduce here a simple but instructive toy model, a single particle following a planar uniform circular motion. By computing probability distributions of the two components of the velocity in this toy case, we will clarify the agreement of the two peaks in the velocity distribution with the ring-like patterns of mosquito motion we observed.
Our toy model is a single particle moving along a circular path, of radius r, at a constant angular speed, ω, in the counterclockwise direction. The particle motion is periodic with period T = 2π/ω. For consistency with the results shown in the main article, we assume that the particle moves on the xz-plane, see Fig.SI1. We measure its angular position, α, as the angle between the x-axis and the radius passing through the particle. We assume that at time t = 0, the particle lies on the x-axis, so that α(t) = ωt, which is linear in time, with α going from 0rad to 2πrad in the time interval [0, T ], see Fig.SI1a.
The 2D position of the particle, ⃗ x(t) ≡ (x(t), z(t)) is described by the following equations of motion: The two components, x(t) and z(t) highlighted in blue and red respectively in Fig.SI1d, are sinusoidal in time, with an amplitude equal to r and a period equal to T . The derivatives in time of Eq.1 give the velocity vector, ⃗ v(t) ≡ (v x (t), v z (t)), of the particle, which is then defined as: As noted for the particle's position, the components of its velocity are also of a sinusoidal form with period T , but with amplitude ωr, see Uniform circular motion. a) A single particle moves, in the xz-plane, along a circular path of radius r at a constant angular velocity ω, corresponding to a revolution period T = 2π/ω. The angular position, α, is measured as the angle between the x-axis and the radius passing through the particle, in the counterclockwise direction. At time t = 0 the particle lies on the x-axis, highlighted with a yellow circle. The red circle represents the position of the particle at a generic time t. In the inset: alhpa is linear in time, going from 0rad to 2πrad during the time period T .b) In the v x v z -plane, the velocity of the particle lies on a circle. The yellow circle represents the velocity of the particle at time t = 0, the red circle its velocity at a generic time t.
In the inset: the angular velocity of the particle is constant. c) The purple line represents the angle α as a function of time, relative to the purple scale on the bottom horizontal axis. The saw tooth shape is due to the equivalence between α and α + 2π. The black line filled in light purple represents the probability distribution of α, relative to the scale on the black horizontal axis at the top. The angle is uniformly distributed, with the probability distribution constant and equal to 1/2π. d) The two components of the position of the particle, as a function of time. The blue line represents the x coordinate, the red line the z coordinate. They are both sinusoidal functions. e) The two components of the velocity of the particle, as a function of time. The blue line represents the x coordinate, the red line the z coordinate. As noted for the position, both components of the velocity also have a sinusoidal trend in time. f) The blue and red lines represent the two velocity components as a function of time, relative to the scale on the bottom horizontal axis. The black line filled in grey represents the probability distribution of the two components of the velocity, which are not constant because of the flat derivative of the sinusoidal function in proximity of ±1.
position of the particle is x = 1 and z = 0, hence its velocity is such that v x = 0 and v z = 1. In the same way, when the angular position of the particle is equal to α, its velocity vector, on the v x v z -plane corresponds to the vector at an angle α + π/2.

Velocity probability distribution
In our toy model, the particle moves at a constant angular speed ω, therefore its angular position α is linear in time, α(t) = ωt, over a period of time T , and displays the periodic saw tooth trend shown with the pink line in Fig.SI1c, over an interval of time longer than T . Its probability distribution, depicted as the light pink solid graph in Fig.SI1c, is constant and equal to 1/2π. On the other hand, the two components of the velocity, v x and v z shown with the red and blue lines respectively in Fig.SI1e, follow a sinusoidal trend in time. Their probability distributions, represented in Fig.SI1f as a solid light grey graph, are equal: they both present a minimum at 0 and two symmetric maxima at ±rω. Therefore, the uniform probability distribution of the angle is not reflected in a uniform distribution of the two velocity components. To better understand how this discrepancy arises, in the following we will first derive the analytic form of the distribution of the velocity components, and then give an intuitive interpretation. Finally, we will show that this distribution is compatible with the mosquitoes' velocity probability distributions presented in the main text. 2/7

Analytic expression
From an analytic point of view, the probability distribution of the angle α, p α (α) shown in Fig.SI2a, is defined as: which is constant and equal to 1/2π in the interval [0, 2π]. The two components of the velocity, defined in eq.2, may be expressed in terms of α as: For the sake of simplicity in the notation, we will first compute the probability distribution of the variable y = cos(α), and only at the end will we generalize the result to give an explicit expression for the two probability distributions of v x and v z , which will differ from p y (y) only by a scale factor. The probability distribution of the variable y, p y (y), and the probability distribution of α, p α (α) are related by the following equation: where dα and dy represent the differential of the two variables α and y respectively, and the factor 2 is needed to take into account that each value of the variable y corresponds to two distinct value of α, see Fig.SI2. Eq.5 is equivalent to: where dα/dy is the derivative of the variable α with respect to y, dy/dα is the derivative of y with respect to α, | · | denotes the absolute value, and the superscript −1 indicates the inverse. So that |dy/dα| −1 = |1/(dy/dα)|.
The derivative of y = cos(α) with respect to the angle α, is given by: Substituting eq.7 in eq.6, we obtain the analytic expression for p y (y) p y (y) is the distribution plotted in Fig.SI2. It has a minimum at 0, where it is equal to 1/π, and it grows symmetrically when y goes to 1 or −1. Note that we find eq.8, also when considering the variable y = sin(α), noting that in this case dy/dα is equal to cos α.
Therefore v x and v z have the same probability distribution that differs from p y (y) in eq.8 only by the overall scale factor rω:

Intuitive interpretation
The analytic expression for the variable y = cos α in eq.8 clearly shows that, a uniform distribution of the angle does not correspond to a uniform distribution of sine and cosine. Here we want to give an intuitive interpretation of this result. We will refer to the cosine function, but the same interpretation holds for the sine function. Why is p y (y) not constant? The cosine, as a function of the angle α, has a flat derivative when its value is close to ±1 (α = 0, π), while it has a derivative close to 1 when its value is close to 0 (α ∼ ±π/2). Therefore, in proximity of α = π/2, the cosine function is almost linear, and uniformly distributed values of the angle will remain almost uniform: an interval of size ∆α in the angle α will correspond to an interval ∆y π/2 in the y variable that is almost equal to ∆α, see Fig. SI2. In contrast, when α is close to 0, the cosine function is not linear because its derivative is almost 0, and the same interval of size ∆α is squeezed into a range, ∆y 0 , that is much smaller than ∆α, hence the density at ∆y 0 is higher than at ∆y π/2 . This higher density corresponds to a higher probability of cos(α) in proximity of ±1. . b) the probability distribution of y = cos(α) is not constant, even when α is uniformly distributed. c) y has flat derivative when it is close to ±1 (α = 0, ±π). The flat derivative implies a low sensitivity of y to a change in α. d) y as a function of α is linear when α is close to ±π/2. An interval of size ∆α, centered on −π/2 is transformed, through the cosine function, to the interval ∆y π/2 , which is almost equal in size to ∆α. Because of the flat derivative of the cosine when α = 0, the same interval of size ∆α centered on 0 is squeezed into the interval ∆y 0 . ∆y 0 has a higher density than ∆y π/2 . e) In orange: the variable y as a function of α (corresponding to the bottom orange horizontal axis). In black filled in light green: the probability distribution of y (corresponding to the top black horizontal axis). The regions where the derivative of y is flat correspond to regions where the density, and hence the probability, is high.

Consistency with real data
In Fig.4 of the main text, we plot the probability distribution of the three components of the velocity at a group level, obtained by aggregating mosquito velocities from 19 single-sex male swarms. In the x and z directions the probability distributions display a particular shape, with a minimum at 0m/s and a symmetric double peak at ±0.5m/s. We find the same double peak distribution in Fig.6, where we show the evolution in time of the z component of the velocity and the corresponding probability distributions for three different males. This trend in the probability distribution is consistent, although affected by noise, with the velocity distribution found in the toy model above, see Fig.SI1e. Experimental data are naturally affected by noise that partly comes from unavoidable 3D reconstruction inaccuracies, but which is mostly biological. Mosquitoes do not perform perfectly circular movements, they do not move at a perfectly constant speed, and they alternate circular movement and straight paths. Moreover, mosquitoes display fluctuations in velocity that are not present in our toy model, where the particle moves at a constant angular speed. All these sources of noise combine, through complicated mechanisms, to transform the shape of the velocity probability distributions, from the one displayed by the toy model to that observed in real swarms. The two sharp edges at the sides of the toy model distribution are replaced by the two distinct peaks in real swarms, and its rigid range of velocity, which cannot exceed ±rω, is replaced by larger range due to mosquito fluctuations of the velocity components. Despite these differences, the toy model is consistent with the data, whose double peak distribution is therefore the effect of the mosquitoes' pseudo-circular motion.

Velocity autocorrelation function
In Fig.6 of the main text, we plot autocorrelation functions of three individual mosquitoes, and use their periodicity to confirm that mosquitoes' motion is characterized by ring-like patterns on the xz-plane. Velocity autocorrelation gives a measure of how much the direction of motion of a particle changes in time, and is a convenient tool for finding periodic and repetitive patterns in the particle path. In order to determine the velocity correlation function, we start by computing c(t 0 ,t) defined as: where · represents the scalar product. c(t 0 ,t) is equal to v 2 when ⃗ v(t 0 ) and ⃗ v(t 0 + t) are perfectly parallel, equal to 0 when the two vectors are orthogonal, and equal to −v 2 when they are anti-parallel, see Fig.SI3. The autocorrelation at time t, C(t), is given by the average over t 0 of the quantities c(t 0 ,t): so that C(t 0 ) = v 2 . In order to compare autocorrelation functions, of objects moving at different speed, we normalized the correlation function by its value in t = 0:Ĉ(t) = C(t)/C(0) in such a way thatĈ(0) = 1.

Fig.SI 3. Toy model -Autocorrelation a) -f)
The purple arrow represents the velocity vector at time t 0 , ⃗ v(t 0 ), the green arrow the velocity vector at time t 0 + t, ⃗ v(t 0 + t). For the sake of simplicity, we assume that they are both unitary vectors. a) ⃗ v(t 0 ) and ⃗ v(t 0 + t) are parallel. Their scalar product is equal to 1, hence c(t 0 ,t) = 1. b) ⃗ v(t 0 ) and ⃗ v(t 0 + t) are orthogonal. Their scalar product is equal to 0, hence c(t 0 ,t) = 0. c) ⃗ v(t 0 ) and ⃗ v(t 0 + t) are anti-parallel. Their scalar product is equal to −1, hence c(t 0 ,t) = −1. d) -f) The purple circle represents the particle at time t 0 , green circle at time t 0 + t. d) At time t = 0, and t equal to multiples of the revolution period T , The autocorrelation function of a particle moving with uniform circular motion,Ĉ(t), is a sinusoidal function, preciselyĈ(t) = cos(ωt), with amplitude 1 (regardless of the radius of the circle along which the particle moves) and with the period equal to the revolution period of the particle.

Real data
When working with experimental data, time becomes a discrete variable defined by the camera frame rate, i.e. the frame rate ( f ps) fixes the minimum time step ∆t = 1/ f ps. Velocity vectors cannot be obtained by computing the derivative of the position vectors, but they are approximated.
From the reconstructed 3D trajectories of individual mosquitoes,⃗ x, we compute individual velocities,⃗ v, from two consecutive frames: where t j = j∆t denotes the instant of time corresponding to the j-th frame, j = 0, · · · , J − 1 with J representing the last frame of the acquisition. Note that we are computing the velocity at frame j from the position at frame j and j + 1, hence we cannot compute the velocity for the very last frame where trajectories are defined.

5/7
The autocorrelation function is defined in discrete time, as for the velocity vectors. The average over time defined in Eq.11 is computed as the following finite sum: with t k = k∆t, J is the last frame of the acquisition, and T max = J∆t is the time corresponding to the last frame. It is worth noting that a value equal to 0 in the autocorrelation, may come from different situations. One possibility is that c(t 0 ,t) = 0 for each t 0 , hence the average over t 0 gives C(t) = 0. This is the case at t = T /4, in the uniform circular motion described above. But there is also a second possibility, more likely to occur with real data: c(t 0 ,t) are not null, but their average is equal to 0. In this latter case, the null value of the C(t) has nothing to do with the orthogonality between ⃗ v(t 0 ) and ⃗ v t 0 +t . C(t) equal to 0 is then the evidence that there is not a typical value for the mutual orientation between ⃗ v(t 0 ) and ⃗ v(t 0 + t), which are then uncorrelated. It is also important to note that for theĈ(t) to be equal to 1 or −1, all the c(t 0 ,t) have to be equal to the normalization factor (all with positive or negative sign). This is a very particular condition, which occurs with the uniform circular motion for t = kπ. But with real data it generally occurs only for t = 0, because of the normalization that actually imposes the conditionĈ(0) = 1.
In Fig.6 of the main text, we show autocorrelation functions of three individual mosquitoes. In order to check whether or not mosquitoes perform ring-like patterns in horizontal planes, we restricted the velocity autocorrelation function to the x and z components of the velocities, in such a way that c(t 0 ,t) is defined as: These autocorrelation functions, restricted to the xz plane, clearly show a typical periodic trend with a period of about 1s. The three plots are quite similar to each other, and are compatible with pseudo-circular motion. The amplitude of C(t) decreases with time, while it would be constant and equal to 1 if the particle were moving in a perfect circle. This discrepancy is due to biological noise, which causes mosquitoes to alternate between circular and straight paths, rather than executing perfect circles, with variations in the radius of the circle and the period of revolution.